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Abstract. We examine the dynamical properties of an exclusion 
process with creation and annihilation of particles in the framework 
of a phenomenological domain-wall theory, by scaling arguments 
and by numerical simulation. We find that the length- and the 
time scale are finite in the maximum current phase for finite 
creation- and annihilation rates as opposed to the algebraically 
decaying correlations of the totally asymmetric simple exclusion 
process (TASEP). Critical exponents of the transition to the TASEP 
are determined. The case where bulk creation- and annihilation 
rates vanish faster than the inverse of the system size N is also 
analyzed. We point out that shock localization is possible even for 
rates proportional toAf~",l<a<2. 
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1. Introduction 

Self-driven many-particle systems have been extensively investigated recent 
years. The ongoing research interest in this kind of systems is both 
conceptual [1-5] and motivated by many important applications in different 
fields [6-13]. Stochastic models of self-driven many particle systems 
have been used in order to describe vehicular transport [6,7], pedestrian 
dynamics [7], intracellular transport [8-10] and many other problems [11- 
13]. The common feature of these models is the steady input of energy, 
which leads to generic non-equilibrium behaviour [13]. This feature implies 
that the standard methods applied in equilibrium statistical mechanics 
have to be generalized in order to handle such systems. First steps in this 
direction have been made: Exact solutions of the stationary state have been 
obtained [14, 15], a quantum formalism was established, which rewrites the 
master equation as a Schrodinger equation in imaginary time [1] , Yang-Lee 
zeroes have been introduced in order to describe non-equilibrium phase 
transition [4] and for a number of transport models a free-energy formalism 
has been developed [3]. 

A question, which naturally arises, is how the properties of such 
systems are influenced by the presence or absence of conserved quantities. 
This issue has been studied by the example of the totally asymmetric simple 
exclusion process (TASEP) [16,17], which is the most simple non-trivial 
driven many particle model. TASEP with particle reservoir which allow 
for particle exchange in the bulk was essentially introduced by Willmann 
et al. as a model of a limit order market [18]. In this paper we study this 
model in the form as defined by Parmeggiani et al. [19] therefore we refer 
to it as PFF model in the following. This model can be viewed in the above 
sense as the grand-canonical counterpart of the TASEP and was motivated 
by the motion of molecular motors, which move along one-dimensional 
filaments [8-11]. It describes correctly the stochastic and biased motion of 
particles, the discrete structure of the filaments as well as the finite length 
of the path between attachment and detachment of a motor, which can be 
tuned by adapting the capacity of the bulk reservoir. If one considers the 
PFF model with open boundaries the particle exchange in the bulk may 
lead to the loc;alization of an interface separating a high- and a low-density 
domain in the bulk. In addition to this, the structine of the phase diagram 
differs strongly from that of the TASEP with particle conservation [19-23]. 
In the present work we investigate the dynamical properties of this model 
by phenomenological and numerical methods. 

The article is organized as follows. In the next section the model is 
re-introduced and the most important features of its stationary state are 
reviewed. In section 3 we examine the dynamical properties of the model in 
the shock phase by means of studying the density-density autocorrelation 
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function which we relate to a phenomenological theory of domain-wall 
motion. The phase diagram of the model is given in the case of vanishing 
total capacity of the bulk reservoir. Phenomenological results are then 
compared to direct simulations of the model. In section 4 we discuss the 
length and time scales in the maximum current phase by means of scaling 
arguments and numerical simulation. A summary and discussion of the 
results follows in the final section. 

2. The model 

The PFF model [18, 19] is defined on a one-dimensional lattice of N sites, 
each of which can either be empty (tj = 0) or occupied by a single particle 
(tj = 1). In the bulk of the system particles interact via asymmetric 
exclusion dynamics, i.e. particle on site i jump to the neighbouring site 
i + I with rate 1 provided it is empty. Boundary sites are coupled to 
particle reservoirs, which realize in- and output rates, whereas bulk sites 
change particles with a bulk reservoir. To be concrete: Particles enter the 
system randomly on site 1 with attempt rate a, and they can leave it on 
the Lth site with rate (3. At bulk sites 1,2,. . . ,L particles can attach with 
attempt rate lua and detach with rate ujd- 

The stationary density (rj) and current profiles (ji) = (rj(l — Ti)) have 
been recently studied by a continuous mean field approximation [19-21]. 
The attachment and detachment rates were taken to be proportional to 
1/N, i.e. LOA = ^a/N, ujd = ^d/N, where ^a and ilo are constants, and 
the thermodynamical limit N oo was considered. The lattice constant 
was rcscalcd according to b = 1/N so the spatial coordinate x = i/N 
becomes a continuous variable in this limit. Neglecting the density-density 
correlations and the spatial derivatives higher than the first-order one, one 
obtains the following equation for the stationary density profile p{x) [20] : 



where K = to a /too- The stationary density profile p{x) can be constructed 
from the flow- field of (1). In order to adapt the solution of (1) to the 
boundary conditions one has to integrate (1) from the left (p(0) = a) and 
right boundary (p(l) = 1 — (3), respectively. This leads to the following 
implicit expressions for the density profile [20] 



(1 - 2p)/ - Qd [K-il + K)p] = 0, 



(1) 



1 K -1 K -{1 + K)p_ 



X 



nD{l+K){l + K) K-{1 + K)a 
2(P- - a) 



+ 



l-x = 



1 K -1 K-{1 + K)I3 



(2) 



nD{l + K){l + K) K-{1 + K)p+ 



Dynamics of an exclusion process.. 



4 



2(/?-p+) 

where (3 — 1 ~ [3 and p_ denotes the solution of (1) obtained from 

integration from the left (right) boundary. The selection of the left or 
right solution is realized by means of characteristics, which determine the 
velocity of discontinuity of the density profile p{x) [20]. If the velocity of 
this so-called shock is finite for any position in the bulk x, it is driven out 
of the system and the stationary density profile is continuous. Contrary to 
this, a localized shock is observed if 

P-{Xa) + P+{Xs) = (3) 

holds for a particular position Q < Xg < 1. The profiles constructed in the 
above way are believed to be exact in the N ^ oo limit [20, 21]. 

In [20] it has also been pointed out that the leading finite-size 
corrections of the density profile in the shock regime can be obtained by 
applying the domain-wall theory for the dynamics of the shock, which was 
originally developed for the TASEP [24] and has been recently generalized 
to models without particle conservation [20, 21]. The idea of this approach 
is to describe the stochastic motion of the domain wall by a random walk 
with hopping rates determined by the particle current in the low- and high- 
density domain. In the case of the TASEP the hopping rates arc constant, 
as the current is constant due to particle conservation. Contrary to this, 
for the PFF model one observes nontrivial current profiles j{x), which lead 
to position-dependent hopping rates: 

p+{x)-p-{x) p+{x)-p- (x) 

where p±{x) and j±{x) = p±{x){l — p±{x)) are the density and the current 
in the high(+) and low(-) density domain, respectively. The potential 
landscape governing the motion of the walker has a minimum at Xg , which 
we will refer to as equilibrium shock position. Previous analysis of the 
TASEP has shown that the random walk picture for the domain wall motion 
gives a correct description of time dependent phenomena as well [26-29]. 
So we believe that this phenomenological description is appropriate also 
for the dynamical properties of the PFF model. 

3. The shock phase 

3.1. The case u)a,d = ^a,d/N 

In this section we discuss the dynamical properties of the model for 
parameter combinations where the density profile has a discontinuity 
separating a high- and a low-density domain. 
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First we consider the case where the total "capacity" of the bulk 
reservoir is comparable with that of the boundary reservoirs, i.e. ujj^ = 
Ua/N, Uj3 = Qjj/N. The density profile of this model was thoroughly 
studied and the parameter regime where the system exhibits a shock is 
known (sec the phase diagram in [20]). 

In order to establish the relevant time scale, we consider the stationary 
(density-density) autocorrelation function C{i, t) = (Ti(O)Ti(t)), whore (. . .) 
denotes the average over the stationary ensemble. One expects that the 
dominant dynamical mode in a finite system, which determines the long 
time behaviour of temporal correlations in the vicinity of the equilibrium 
shock position Xg , is the stochastic motion of the domain wall as opposed 
to "microscopic" processes. Thus, for large N the local density at a given 
time t is appropriately described by the function 

T{x,t) = p+{x) + ip-{x) - p+{x)) 9{m - , (5) 

where 6{x) is the Heaviside function and ^{t) is a random walk with steps of 
length 1/N and with hopping rates given in (4). The potential well which 
the walker is trapped in is well approximated by a harmonic potential in 
the vicinity of its minimum Xg- Considering the continuous description of 
the random walk the Fokker-Planck equation reads [32] 

dP V d , D d^P 

where y = x — Xg is the deviation from the equilibrium shock position. The 
constants V and D, which characterize the shape of the potential, are given 
by 

^ _ dwi{Xs) _ dWrjXs) 

dx dx ' 

D = Wl{Xs) — Wr{Xs). (7) 

For the PFF model we have V = ^Ia + ^d- Equation (6) is the Fokker- 
Planck equation of the Ornstein-Uhlenbeck process [32] , and has the time- 
dependent solution: 

P{y,t) = P{y,t\yo,Q) 



2itD 



-(1-e 



-1/2 



X exp 



N'^{uja+u}d) 

N^joJA + ujp) {y - t/oe-('^^+"°)*)^ 
2D 1 — e-2(wyi-i-WD)t 



.(8) 



For large system sizes the harmonic approximation is expected to give an 
appropriate description of the shock dynamics, since the localization length 
of the walker increases only sub-extensively (~ N^^^). 
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In the framework of the above phenomenological picture the 
autocorrelation function is given by 

C(t/,t) = (r(y,0)r(y,t)) 

= p-{y)[2p{y) - p.{y)] + [p+{y) - p.{y)]'l{y,t), (9) 

where 

i{y,t)^{em)-ynm-y)) 

/y rv 
dXo dxP{x + Xs,t\xo + Xs,0)Pstixo + Xs). (10) 

Inserting now (8) into (10) wc obtain 



-{x+yf 



xerfc I (^(^ + v) - y)j dx, (11) 

where T = Q-(.'^A+<^n)t ^^^^^ erfc(a;) is the complementary error function. 

The integral I{y,t) can be evaluated only at the equilibrium position 
of the shock, i.e. for y = 0, where we obtain the well-known result 

7(0, t) = 1 + -L arcsin[e-('"'^+'^-°)*]. (12) 

For y ^ I{y,t) cannot be calculated analytically, however, it can be 
expanded for short ((w/v + L^'_D)t <C 1) and long times {loa +cOD)t '> 1). In 
the latter case the asymptotic form of I{y, t) is given by: 



i{y,t) 



1 




l + erf 



N'^{uja+ujd) 



2D 



+ 0{e 



-2{u)A+lt}D)t 



) ■ (13) 



This expression shows that the motion of the domain wall introduces a 



time scale t = 



AT 



which is proportional to the system size. 



Note that this time scale r is independent of the position y. However, the 
domain wall contribution to the true correlation function is relevant only 
in the region \y\ N~^^'^ and its amplitude is exponentially suppressed 
when leaving the equilibrium shock position. 

It is also interesting to discuss the short time behaviour of I{y,t), i.e 
the case {wa + 0JD)t <^ 1. Here, the expansion of (11) yields 



i{y,t) 




N^{uJA+iOD) 



2D 



— e 



JvVA±iiDl,,2 1 /wa + OJd 

2D y . ' 



+ , 



(14) 



Dynamics of an exclusion process.. 



7 



The leading order correction is thus proportional to t^/^, similarly 
to the TASEP with parallel dynamics, where the space- and time- 
dependent correlation function is exactly known [25] . This asymptotic form 
corresponds to a /~^/^ power spectrum of the local density fluctuations as 
it was found in the case of the TASEP [26]. 

3.2. The case of vanishing total capacity of the bulk reservoir 

We now turn to discuss the case where to a and lod vanish faster than 1/A''. 
In [19, 21] it was then argued that the effect of bulk reservoir is negligible 
in the thermodynamic limit, and the system will behave as the TASEP. 
We have found that this scenario is correct for any parameter combination 
except the line a — P < 1/2. Consider at this particular line the general case 
where attach and detach rates scale as iJA,D = ^a,d/N". We claim that 
a shock in the density profile is still observed in the thermodynamic limit 
whenever 1 < a < 2. Substituting ua.d = ^a,d /N"" into the calcidations of 
the previous subsection one obtains that the width of the shock grows with 
the system size as ^ = NJS.X r-^ Ni , whereas the time scale as r ~ N°'. 
Therefore, as long as 1 < a < 2, the localization length of the domain 
wall increases only sub-extensively, and an unbounded motion (as for the 
TASEP with a = jS < 1/2), is only observed if a > 2. This is also true 
for the related time scale r. The time scale in the TASEP is known to 
diverge proportionally to N'^ [26,27,29]. Therefore if 1 < a < 2 the time 
scale is determined by the attach and detach processes (r ~ N"), while 
TASEP-like behaviour (r ^ N"^) is observc;d only if a > 2. 

In the case 1 < a < 2 the average density takes the values a in the 
low-density domain and 1 — a in the high-density domain if A'' — > oo. Thus, 
the location of the discontinuity in the N ^ oo limit can be obtained by 
substituting p-{x) = const = a and p+{x) = const = 1 — a into (1). 
Equation (3) then yields 

_ l-a{l + K) 
'^^-(l+i^)(l-2a)- ^^^^ 
Note that Xs is independent of a. Thus, a shock can be found in the system, 
whenever 

a. 1 — ex 
<K < . 



1 - a 

If > K (K > ^^^) the system is in the (low-) high-density phase. 

For a ^ j3 and a > 1 it is not possible to balance the current in the 
low- and high-density domain, i.e. to fulfill relation (3) for any < a; < 1 
in the N oo limit. Hence the shock is driven out of the system and the 
density profile of the TASEP is recovered. 

The description of the domain- wall dynamics by a random walk implies 
that the dynamical exponent (defined as r ^^) is 2; = 2. The same 
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exponent can be observed for the TASEP at the transition hne a = < 1/2, 
where the length- and time scale are infinite in the limit N ^ oo, whereas 
for finite systems they scale as ^ ~ iV and t ~ A^^, i.e. z = 2. In the case of 
the TASEP this divergence is restricted to the transition line a = (3 < 1/2, 
because the domain wall motion is biased if a ^ 0. In the case of the PFF 
model the position dependent hopping rates introduce a localisation length 
which grows 

sub-cxtcnsivcly in the whole shock phase. Here, the transition 
between shock- and low-density phase or shock- and high-density phase 
is not a localization-delocalization transition, but it simply means that 
the equilibrium shock position moves to the system boundaries. At the 
transition line the random walker is trapped in a potential with hard- 
core repulsion for y < (HD-S transition) or ?/ > (LD-S transition) 
respectively. The complementary part of the potential landscape can still 
be described by the harmonic approximation. 

Finally we mention that the domain wall contribution of the connected 
autocorrclationfimctionCconn(y, t, A'') = (r(y, 0)r(y, t)) — (t(?/,0))^ has the 
following scaling form: 

Cconn{y,t,N) = C{y^N^-'',tN-''). 

We stress that the contribution to large time behaviour of the 
autocorrelation function, induced by the domain-wall movement is relevant 
in a finite system only at sites |y| <^ N^~^. 

3.3. Numerical results 

In this subsection we compare the phenomenological results obtained in the 
previous parts of this section with results of direct numerical simulation of 
the model. In the numerical investigations we have mainly considered the 
case a = 1, = Q.d = where the equilibrium position of the shock Xs 
and the diffusion constant D can be explicitly given. Equations (3) can be 
solved for p+(x) and P-{x) and condition (3) yields Xs = + \. The 
diffusion constant is given by D = i(l/A — A)) where A = 1 — a — (3 — Q, 
is the height of the shock [20] . 

The connected autocorrelation function was computed for system sizes 
A'^ = 201 — 2001. Results at the equilibrium position of the shock [y = 0) 
are shown in figure 1 and figure 2. As can be seen in the figures the 
phenomenological predictions arc in a good agreement with the numerical 
results except for short times, and the accuracy of the phenomenological 
description improves for larger system sizes, as it was found in the case 
of the stationary density profile [20]. A rapid decay of the autocorrelation 
function can be observed on a time scale t ~ 0(1), which is connected to the 
"self-correlation" of particles, i.e. the contribution of particles which are 
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Figure 1. The connected autocorrelation function measured at 
y = for different system sizes and with parameters a = 0.1, /3 = 0.1 
and Q = 0.1. The solid line is the phenomenological prediction given 
by (9) and (12). 



not updated during time t. This "microscopic" time scale is related to the 
particle current and can be identified as the finite time scale in the high- and 
low density domains of the TASEP [27,29]. The above phenomenological 
picture apparently does not account for this contribution, nevertheless the 
crossover to the regime, which is dominated by the motion of the domain 
wall, takes place at rather short times of 0{\/N). 

In figure 2 results for sites away from the equilibrium shock position 
{y 7^ 0) ^re presented. The accuracy of the phenomenological theory is 
less satisfying as leaving the equilibrium shock position. This discrepancy, 
which increases with larger distances from the equilibrium shock position, 
may be related to the anharmonicity of the potential-well. 

4. Maximum current phase 

The existence of a maximum current phase in the presence of bulk particle 
exchange is restricted to the case oja = <jJd = w (sec [20]). Furthermore, 
we assume that the rate lu docs not scale with the system size, i.e. a = 0. 
Compared to the maximal current phase of the TASEP we expect that 
temporal correlations are reduced, as the typical lifetime of particles is 
finite, and this may introduce a finite time scale. This was also observed in 
[18] where the correlation function and fiuctuations of price increments were 
investigated and the corresponding time scale was found to be Here, 
we study the connected density- density autocorrelation function at the site 
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t/N 



12 



Figure 2. The connected autocorrelation function measured at 
different positions in the chain for system size N = 1001 and with 
parameters as in figure 1. The phenomenological curve for y ^0 was 
obtained by integrating (11) numerically. 



in the middle of the chain. Since the vaUdity of the phenomenological 
treatment presented above is restricted to the shock phase, we resort 
to scaling arguments and numerical simulations. To our knowledge this 
correlation function has not yet been investigated in the maximum current 
phase even in the case of particle conservation in the bulk. Therefore we 
examine the dynamical correlations in the TASEP first. 

Let us assume that the connected autocorrelation function is a 
homogeneous function of its variables N and t. Rescaling then the lengths 
by a factor b it transforms as 

C{N, t) = b-''C{N/b, t/b^/^), w = (16) 

where we have used that the dynamical exponent is ^; = 3/2 in the 
maximum current phase [29,30]. The scaling dimension x can be guessed 
as follows. The fluctuations of the total particle number scale with the 
system size as ~ N^^"^ [33,34]. As a consequence the fluctuations of the 
local density scale as ^ N~'^/'^. Since the autocorrelation function contains 
a product of two local density operators, we have x = 1. Setting now 
t = 6^/2 in (16) and taking the limit iV ^ 00 we obtain C{t) ~ 
whereas choosing b = N wc get 

C{N,t) = N-^<^>{t/N^/^), (17) 

where the scaling function $(a;) behaves as 

i>(a;) ~ x~'^/^ X < 1 

^{x) ~ a; > 1. (18) 
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2 4 6 8 



ln(t) 

Figure 3. The connected autocorrelation function in the maximum 
current phase (a = /3 = 0.5) of the TASEP {uja = ujo = 0) for 
different system sizes. 

Numerical results for the autocorrelation function are shown in figure 3, 
indicating an algebraic decay with an exponent compatible with 2/3, and 
the scaling plot in figure 4 is in accordance with (17). 

For finite uj and infinite system size one expects the following behaviour 
of the autocorrelation function when time is rescaled by a factor b: 

C{t,uj) = b-'^/^C{t/b,Ljb), 1/N = 0. (19) 

Setting b = we have 

C{t, oj) = oj^^^^itcu), 1/N = 0, (20) 

where the scaling function ^{x) behaves as $(x). Numerical results for 

finite Lu are shown in figure 5 and figure 6. These arc in agreement with 
(20), showing a cut-off at a time scale r ~ w"^ which is the typical lifetime 
of particles [18]. 

Next, we determine the length scale ^ in the infinite system. In 
order to establish ^ we have studied the asymptotical decay of the density 
profile. For the TASEP the density profile is known to approach its bulk 
value 1/2 algebraically [15, 17]. Contrary to this, the particle creation and 
annihilation processes lead to an exponentially fast asymptotical approach 
to the bulk density. To see this, consider the second-order mean field 
equation for the stationary density profile [19] 
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Figure 4. The scaled connected autocorrelation function in the 
maximum current phase (a = j3 = 0.5) of the TASEP (loa = ^^d = 
0). 
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Figure 5. The coinicctcd autocorrelation function in the maximum 
current phase (a = j3 = 0.5) for different rates u). The size of the 
system is AT = 257. 
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Figure 6. Scaling plot of the connected autocorrelation function in 
the maximum current phase {a = j3 = 0.5) according to (20). 



The solution of (21) in the asymptotic region where p{x) — ^ <C 1, has the 
form p{x) - 5 ~ 6"^^/^-/. Putting this into (21) yields = The 
scale of the decay S^mf can then be identified as a finite length scale of the 
system. However, the true behaviour of ^ is not recovered by the mean 
field approximation as we show in the following. 

Considering w as a control parameter, the point lu = (TASEP) can 
be regarded as a critical point of the PFF model, where length- and time 
scale diverge. We expect that in the vicinity of the critical point, i.e. for 
the PFF model with w <C 1 the dynamical exponent is the same as strictly 
at criticality, i.e. that t ~ ^^^^ holds. Comparing this with r ~ uj~^ we 
obtain 

^ ~ uj-", V = 2/3. (22) 

In this context the exponent ly plays the role of the correlation length 
exponent. The correctness of (22) is supported by numerical results for the 
density profile shown in figure 7 and figure 8. 

Thus, we have seen that the length- and the time scales remain finite 
in the thermodynamic limit as long as oj does not vanish in that limit, i.e. 
a = 0. If o > 0, apparently both the length- and the time scale diverge 
in the ^ oc limit, as in the case of the TASEP. However, for a < | 
they scale with exponents that are different from those of the TASEP, viz. 
^ ~ A^ 3^, r ~ A^'^, where c = min{a, |}. 
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Figure 7. Density profile in tlie maximum current pliase (a = /3 = 
1) for diflferent rates ui. The size of the system is N=512. 
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Figure 8. Scaling plot of the density profile in the maximum current 
phase (a = /? = 1). The size of the system is N=512. 

5. Discussion 

The analysis of the dynamical properties completes the description of the 
stationary state of the PFF model. In the framework of a phenomenological 
domain-wall theory we could establish a time scale in the vicinity of the 
equilibrium domain wall position, which grows with the system size as 
~ A^". The localisation length of the domain- wall scales as ~ Aft, leading 



Dynamics of an exclusion process.. 



15 



to the dynamical exponent z = 2, which is identical to that of the TASEP 
at the coexistence line a = /3. This is in both cases the consequence of the 
diffusive nature of the dominant dynamical mode. However, in the PFF 
model it is relevant only in a subextensively growing region whereas in the 
latter one in the whole system. Aside from the shock region the time scale 
is finite, and is related to the inverse of the particle current. Contrary to 
the TASEP the current is position dependent, a feature of the model which 
is reflected by the position dependence of the time scale. This position- 
dependent microscopic time scale is observed in the high- and low-density 
phases, as well. We note that the transition from the shock- to high-density 
(or low-density) phase manifests itself simply in the equilibrium domain- 
wall position leaving the system and it is not a delocalisation transition 
as in the case of the TASEP. We have pointed out that for a = /? < 1/2 
even a vanishing total capacity of the bulk reservoir (1 < a < 2) is able 
to localise the shock. We mention that the phenomenological calculations 
presented in this work can be carried out for arbitrary systems possessing 
a localised fluctuating domain wall, such as the model recently introduced 
by Rakos et al. [35]. 

In the maximum current phase of the TASEP we have found the 
algebraic decay of the autocorrelation function. By scaling arguments we 
could determine the decay exponent which was found to be 2/3. The 
introduction of noise through bulk particle exchange with a finite rate 
destroys the power-law correlations and the resulting phase is characterised 
by finite length- and time scales. Thus, the TASEP can be regarded as a 
critical point of the PFF model. By scaling arguments we have determined 
the critical exponents which are in accordance with the results of numerical 
simulations. 
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